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ABSTRACT 

There is mounting observational evidence that most galactic nuclei host both supermassive 
black holes (SMBHs) and young populations of stars. With an abundance of massive stars, 
core-collapse supernovae are expected in SMBH spheres of influence. We develop a novel 
numerical method, based on the Kompaneets approximation, to trace supernova remnant 
(SNR) evolution in these hostile environments, where radial gas gradients and SMBH tides are 
present. We trace the adiabatic evolution of the SNR shock until 50% of the remnant is either 
in the radiative phase or is slowed down below the SMBH Keplerian velocity and is sheared 
apart. In this way, we obtain shapes and lifetimes of SNRs as a function of the explosion 
distance from the SMBH, the gas density profile and the SMBH mass. As an application, we 
focus here exclusively on quiescent SMBHs, because their light may not hamper detections 
of SNRs and because we can take advantage of the unsurpassed detailed observations of our 
Galactic Centre. Assuming that properties such as gas and stellar content scale appropriately 
with the SMBH mass, we study SNR evolution around other quiescent SMBHs. We And that, 
for SMBH masses over ~ 10^ M 0 , tidal disruption of SNRs can occur at less than 10"* yr, 
leading to a shortened X-ray emitting adiabatic phase, and to no radiative phase. On the other 
hand, only modest disruption is expected in our Galactic Centre for SNRs in their X-ray stage. 
This is in accordance with estimates of the lifetime of the Sgr A East SNR, which leads us to 
expect one supernova per 10'^ yr in the sphere of influence of Sgr A*. 

Key words: accretion, accretion discs — black hole physics — hydrodynamics — shock 
waves — ISM: supernova remnants — galaxies: nuclei 


1 INTRODUCTION 


There is compelling evidence for a supermassive black hole 
(SMBH) with a mass of 4.3 x 10® Mq in the nucleus of the Milky 
Way, associated with the Sgr A* radio source. The strongest evi¬ 
dence comes from the analysis of orbits of the so-called ‘S-stars’ 
very near this compact object, such as that of t he star S2 with a pe¬ 
riod of only 16 yr and pericentre of ~ 10^ au jSchodel et al.l 2002 , 
200^ Ghez_et^ 20031 : lE^senhauer et aflbOO? Ghez et alj|20oi : 
Gillessen et al.ll2009h. 

Most other massive galaxies contain SMBHs_jMarl£au_et_^ 


20131) , some with masses as high as 10^® Mq iMcConnell et al 


201 ih . The observed fraction of active nuclei is no m ore than a 
few per cent at low redshifts dSchawinski et al.ll2010h . and most 
alactic nuclei house very sub-Eddington SMBHs, like Sgr A* 


Melia & Falckel200l] : lAlexandeil200^lGenzel et al.l2010h . These 
SMBHs are believed to be surrounded by radiatively inefficient ac¬ 
cretion flows (RIAFs), where only a small fraction of the accre- 
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tion energy is carried away by radiation ( llchimariJl977l : [Rees et al.l 
ll982l:lNaravan & Yill994l) . 


In addition to the ubiquity of SMBHs, young stellar 
populations and appreciable star formation rates are corn- 


mon in many quiescent galactic nuclei JSarzi et al. 

2005 

Walcher et al.ll200d:ISchruba et al.ll201ll:lKennicutt & Evans! 

2012 

Neumaver & WalchehoiTI) ^hhis is seen most clearly in the abun- 


dance of early-t yp e stars in the ce ntral parsec of the Milky Way (see 
IDo et al.ll2013a| [a: fLu et alJl2013L for some recent reviews). More¬ 
over, it appears that star formation in the Galactic Centre region 
has been a persistent process that has increased over the past 10® yr 
jFiger et al.l20o3:lFigej26o^ : |Pfuhl et alj201 ih . Over that time, an 
estimat ed > 3 x 10® Mp) of stars have form ed within 2.5 pc of the 
SMBH ( iBlum et al.ll20^:|Pfuhl et aDl201 ll) . 

Continuous star formation in galactic nuclei will regularly re¬ 
plenish the supply of massive stars in these regions. This natu¬ 
rally leads to the expectation of frequent core-collapse supemovae 

^ Evidence for recent star form ation has also been seen around active 
galactic nuclei (AGN; for example. iDavies et al.l200'7h . However, active nu¬ 
clei are not the subject of this study. 



















































2 Rimoldi et al. 


in such environments. As an example, IZubovas et al.l l l2013ll show 
that, per 10® Mq of stellar mass formed in the Galactic Centre, 
approximately one supernova per 10"^ yr is expected for the past 
10® yr. 

Only one supernova remnant (SNR) candidate has been iden¬ 
tified close to the SMBH sphere of influence (SOI): an elongated 
shell known as Sgr A East, at the end of its adiabatic phase. It 
has an estimated age of about 10"^ yr and appears t o be engulf¬ 
ing Sgr A* with a mean radiu s of approximate l y 5 pc l Maeda et al.l 
l2()02l:]Hermstein & Hdl200^: iLee et alJ[200^ : iTsuboi et al.ll2009tl . 
In addition, there are a couple of observations that indirectly point 
towards supernovae in the SOI. The first is CXOGC J174545.5- 
285829 (‘The Cannonball’), suspected to be a runaway neutron 
star associated with the same supernova explosion a s Sgr A East 
dPark et al.ll2005l : lNvnka et aDl2013l : IZhao et ^l2013h . The second 
is the recently discovered magnetar SGR J1745-2900, estimated to 
be within 2 pc o f Sgr A* jPegenaar et al.l2Ql'^ : lKennea et al.l2ol3l : 
iRea et alJl2Q13[) . 


Any supernova exploding in the SOI of a quiescent SMBH 
will expand into a gaseous environment constituted mainly by the 
SMBH accretion flow, whose gas is supplied by the winds from 
massive stars. The density distribution within the flow is there¬ 
fore set by both the number and distribution of young stars and 
the hydrodynamical properties of a radiatively inefficient accretion 
regime. This interplay gives an overall density distribution that is 
a broken power law, for which the break occurs where the number 
density of stellar wind sources drops off. For the Galactic Centre, 
this corresponds to ~ 0.4 pc (for example, l0uataerll|2004l) . 

In such environments, we expect SNRs to evolve differently 
from those in the typically flat interstellar medium, away from the 
SMBH. The density gradients have the potential to distort SNRs 
and decelerate them significantly. Once the expansion velocity falls 
below the SMBH velocity field, the remnant will be tidally sheared 
and eventually torn apart. This can substantially shorten an SNR 
lifetime compared to that in a constant-density interstellar environ¬ 
ment. In turn, this can reduce the expected number of observed 
SNRs in galactic nuclei. 

Since quiescent accretion flows are fed by stellar winds, which 
can be also partially recycled to form new stars together with the 
gas released by supernova explosions, the scenario we consider is 
of a self-regulating environment, where young stars and gas (or, 
in other words, star formation and accretion on to the SMBH) are 
intimately related. This holds until a violent event—for example, 
a merger—drives abundant stars and gas from larger scales to the 
galactic nucleus. Observations and modelling of our Galactic Cen¬ 
tre support this picture. In particular, winds from massive stars are 
sufficient to account for the observed accretion luminosity and ex¬ 
ternal gas feeding is not required (e.g. lOuataerll2()04l : ICuadra et al.l 
l2006h or observed. Furthermore, ther e is strong evidence for the 
recent star formation occurring in situ dPaumard et al.ll20o3) . 

In this paper, we determine the morphology and X-ray life¬ 
times of SNRs, which, in turn, can be used to constrain the envi¬ 
ronment of SMBHs. We develop a numerical method to trace SNR 
evolution and determine their X-ray lifetime. The influence of the 
SMBH on SNRs will be considered first indirectly, through its in¬ 
fluence on the gaseous environment, and then directly, through its 
tidal shear of the ejecta. 

The paper is organized as follows. Section introduces the 
gaseous environments found around quiescent SMBHs. Section |3 
uses analytic methods to qualitatively trace SNR evolution. Sec¬ 
tion |4] describes our numerical method, which allows us to follow 
the evolution of an SNR in an arbitrary axially symmetric gas dis¬ 


tribution. We then specialize it to a quiescent SMBH environment. 
Section|5]outlines the galactic models used for the environments of 
the supernova simulations. Section|^presents our results for SNR 
shapes and lifetimes. Our concluding remarks are found in Sec- 
tion|7] 


2 GASEOUS ENVIRONMENTS OF QUIESCENT 
GALACTIC NUCLEI 


In this section, we outline the expected gas distributions near the 
SMBH in quiescent galactic nuclei. These gas distributions will be 
used as the environment for the SNR model exposited in Sections 
[3] and |4] We will then proceed to scale the general environment 
discussed here for the Galactic Centre to other SMBHs in Sectionj^ 

Quiescent SMBHs are surrounded by RIAFs, which are the 
environments in which the SNR will evolve. RIAFs are relatively 
thick, for which the scale-height, H, is comparable to the radial 
distance, R, from the SMBH (H/R « 1). The mechanisms of en¬ 
ergy transport within the flow vary depending on the model, and 
these variations affect the power-law gradient in density near the 
SMBH. Advection-dominated accretion flow (ADAF) models as¬ 
sume that much of the energy is contained in the ionic compo¬ 
nent of a two-temperature plasma. As the ions are much less effi¬ 
cient radiators than electrons, energy is advected i nto the SMBH 
by the ions before it can be lost via radiation jNaravan et al.l 
1 19951 : iNaravan & Yil Il995ll . Additionally, convection-dominated 
accretion flow (CDAF) models rely on th e transport of energy 
outward via convectiv e motions in the gas JOuataert & GruzinovI 
1 20001 : iBalletaljgIM ). Finally, the adiabatic inflow-outflow rn odel 
(ADIOS: fBiandford & Begelmanlll999l.[2004fBegelmaij|2012ll ac¬ 
counts for winds from the flow that expel hot gas before it is ac¬ 
creted. 

For the region near the SMBH, predicted exponents, ojin, of 
the power law in gas density, p, lie in the range of ojin = 1/2 
to 3/2. The lower and upper limits of u;in are derived from the 
predictions of the CDAF/ADIOS and ADAF models, respectively. 
A drop-off in stellar number density at a radius R = Rh from the 
SMBH would cause a break in the mass density, p, at the same 
radius, since it is the winds from these stars that feed the accretion 
flow. 


The best example of a RIAF is that surrounding Sgr A*. It 
has been extensively studied theoretically and observationally and 
will constitute our prototype. A density distribution from the one¬ 
dimensional analytic model of wind sources has approximately a 
broken power-law sha pe with cJin = 1 inside the density break 
and tJout = 3 outside dOuataertll2004h . Simul ations of stellar win d 
accretion show comparable density profiles iCuadra et alj|2005l . 
Furthermore, the value of Win = 1 is consistent with GRM HD ac¬ 
cretion simulations (for example. iMcKinnev et alj[2012h . Recent 
observations using long integrations in X-ray suggest that a gradi¬ 
ent of Win ~ 1/2 may provid e a better fit to the inner accretion flow 
of Sgr A* l lWang et al.ll2013h . 

We can therefore, generally describe the ambient medium of a 
quiescent SOI with a broken power law for the density of the form: 


p{R) = 


' 

( R\ 

Po 

K,Rh) 


( R\ 

ph 



R ^ Rh 


R > Rh, 


( 1 ) 


for Win G {1/2, 1, 3/2}, Wout = 3, using a reference point for the 
density at R = Ro away from the SMBH. 
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The strongest observational constraint on the density around 
Sgr A* is given by Chandra X-ray measurements at the scale of 
the Bondi radius (_Ro ~ 0.04 pc) of np ~ 130 cm~® (po ~ 
2.2 X 10“^^ g cm~^-. [Baganoff et aljl2003h . The accretion rate 
closer to the SMBH can be further constrained by Faraday rota¬ 
tion m easurements, though the relative error is large jMarrone et al.l 
l2007h . Indeed, we find that fixing the density at 0.04 pc and vary¬ 
ing Win between 1/2 and 3/2 produces a range of densities at small 
radii that fall within the uncertainty in the density inferred from 
Faraday rotation. The radius for the break in stellar number density 
and gas density in the Milky Way is taken to be i?b = 0.4 pc. 


3 EVOLUTION OF REMNANTS AROUND QUIESCENT 
BLACK HOLES: ANALYTIC FOUNDATIONS 

Here, we outline the physics describing the early stages of SNR 
evolution that are of interest in this work. The theory described in 
this section will be used as the foundation of a general numerical 
method to solve the problem, outlined in Section [T] At this point, 
we do not directly take into account the gravitational force of the 
SMBH, but instead just the gaseous environment. The gravity of 
the SMBH can be ignored when the expansion velocity of the SNR 
is much larger than the Keplerian velocity around the SMBH. For 
example, around Sgr A*, at a velocity of 10^ kms“^ gravity can 
be ignored for radii larger than ~ 10~^ pc. The gravitational field 
of the SMBH will be accounted for later, when we consider tidal 
effects on the expanding remnant, which are important only once 
the remnant has slowed down significantly. 

A supernova explosion drives a strong shock into the sur¬ 
rounding gas at approximately the radial velocity of the ejected de¬ 
bris. Typically, it is assumed that a significant amount of the ejecta 
is contained within a shell just behind the shock front (for example, 
IKoo & McKed[T99^ . As it expands, the shock sweeps up further 
mass from the surrounding medium. By momentum conservation, 
the combined mass of the fraction of ejecta behind the shock front 
(Mej) plus the swept-up gas (Ms) must decelerate. The decelera¬ 
tion is considered to be appreciable when the swept-up mass be¬ 
comes comparable to that of the debris, and therefore this ejecta- 
dominated phase holds for Ms <C Mej. 

The subsequent adiabatic expansion of the shock front is mod¬ 
elled with the assumption that losses of energy internal to the rem¬ 
nant are negligible. For this decelerating regime, the Rankine- 
Hugoniot strong-shock jump conditions can yield exact similar¬ 
ity (length scale-independent) solutions for the kinematics of the 
shock front. The evolut ion is determined by its energy, E, and 
the ambient density, p iMcKee & Truelovml995ll . In all of this 
work, we use a canonical value of 10®^ erg for the explosion en¬ 
ergy. In a uniform ambient medium, the adiabatic stage is classi¬ 
cally modelled using the spherica lly symmetric Sedov-Taylor so¬ 
lution (iTavloill 19501 : ISedo\lll959l) . This has self-similar forms for 
the spherical radius and speed of the SNR of R' oc {E/pf^^ 
and V oc {Ejpf^^ respectively, where R' is measured from 

the explosion site. _ 

_ Following the initial work by Sedov and Taylor, iKompaneetsI 

( Il960h developed a non-linear equation from the jump conditions 
that allows self-similar solutions for the shock front evolution 
in certain density stratifications. The original work by Kompa- 
neets considered an atmosphere with exponential stratification, but 
many other solutions have since been obtai ned (see the review by 
iBisnovatvi-Kogan & Silichiri995L as well as [Bannikova et alj|20ld 
and the references therein). Of particular relevance to the gas dis¬ 


tributions in galactic nuclei. iKorvcansk^ i ll9921) —hereafter, K92— 
showed that, with a specific coordinate transformation, a circular 
solution to the Kompaneets equation can be obtained for explo¬ 
sions offset from the origin of a power-law density profile, 

(for u! 2). 

The early ejecta-dominated and late adiabatic stages are well 
characterized by the purely analytic solutions for each stage. In 
between, the solution asymptotically transitions between these 
two limits (this is known as ‘intermediate-asymptotic’ behaviour; 
iTruelove & McKej|l999l) FI The late evolution of the remnant, the 
radiative stage, occurs when the temperature behind the shock 
drops to the point at which there is an appreciable number of bound 
electrons. Consequently, line cooling becomes effective, the ra¬ 
diative loss of energy is no longer negligible, and the speed of 
the shock will drop at a faster rate. For SNRs in a constant den¬ 
sity of n ~ 1 cm~^, the radia tive phase begins at approximately 
3x10^ yr telondin et al.ll99^ . We do not model the remnant dur¬ 
ing this phase, but we will estimate the onset of the transition to the 
radiative stage. 

In the present work, we model SNRs over the first two (ejecta- 
dominated and adiabatically expanding) stages of evolution in a 
range of galactic nuclear environments. The evolution begins with 
a spherically expanding shock, and therefore we do not consider 
any intrinsic asymmetries in the supernova explosion itself. Collec¬ 
tively, any possible intrinsic asymmetries in SNRs are not expected 
to be in a preferential direction, and so they should not bias the 
generalized results presented here. 

The overall geometry of this analysis is laid out in Fig. [T] 
which indicates the main coordinates, distance scales and density 
distributions. The explosion point is at a distance R = a, measured 
from the SMBH (the origin of our coordinate system). The shock 
front extends to radial distances R', measured from the explosion 
point. Each point along the shock is at an angle tp, measured from 
the axis of symmetry about the explosion point. The initial angle 
made with the axis of symmetry of each point on the shock, at 
f —>■ 0, is denoted tpo. 

3.1 End of the ejecta-dominated stage 

In order to estimate where the shock front kinematics appreciably 
deviate from the ejecta-dominated solution, we integrate the back¬ 
ground density field along spherical volume elements swept out 
by the expanding remnant. This provides an estimate of the mass 
swept up from the environment, Ms. The ejecta-dominated solution 
is taken to end when Ms is equal to some specified portion of the 
ejecta mass, Mej. We use a canonical value of IM© for this frac¬ 
tion of ejecta mass. The distance from the explosion point (along 
the coordinate R') at which this occurs is denoted the ‘decelera¬ 
tion length’, L, here (it also known as the ‘Sedov Length’ in the 
standard treatment of SNRs in a uniform p). 

Since our density profiles are not uniform, different directions 
of expanding ejecta will sweep up mass at different rates. In gen¬ 
eral, we must consider a solution for L that depends on ipo, the ini¬ 
tial angle of each surface element of the shock with respect to the 
axis of symmetry (see Fig.[T]l. We therefore determine the value of 
L{tpo) corresponding to small surface elements of the shock front. 
When the explosion occurs close to the SMBH, the solution is ex¬ 
pected to converge to that of an integral over a sphere, due to the 

^ For an illustration of this transition, see fig. 2 of iTmelove & McKed 
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shock segment 



Figure 1. Basic geometry of the problem. The supernova occurs at a point S, a distance R = a away from the SMBH, which is located at the origin, O. The 
shock front extends to distances measured radially from the explosion point S by the coordinate R'. The angle made by a point on the shock, measured from 
the 0 = 0 axis about the explosion point, is denoted ip. Each point on the shock has an initial angle ip{t 0) = ipo- The entire density distribution p{R) can 
be characterized by: the choice of the inner gradient arm (defining the density within the shaded circle), the outer gradient cJout, the reference density po (at a 
reference radius Ro), and a break at between the gradients tUin and tUout- 


spherical symmetry of the background density^ Therefore, as a ref¬ 
erence, we also find the radius L of the sphere whose volume en¬ 
closes Ms « Mej . 

The explosion occurs at a distance R = a from the origin. For 
a single power-law stratification, we use the explosion point for the 
reference density, po = p(a) = Pa, such that 

P{R) = Pa . (2) 

We consider a small surface element of the SNR at an angle ipo 
over an infinitesimal solid angle. In a single power-law stratification 
with the form of equation the length L (ipo) can be estimated 
from the mass integrated through R' at a given angle ipo : 

= Mej. (3) 

Note that we are integrating over the coordinate R' that extends ra¬ 
dially from the explosion point, but that the density varies radially 
with the coordinate R as measured from the SMBH. For integrals 
over a broken power-law density, the density break adds complica¬ 
tions to the integrals analogous to equation ID- The solutions are 
discussed further in Appendix lAl 

These methods for estimating the deceleration length provide 
a means for testing the level of asymmetry and distance scales in the 
ejecta-dominated stage of evolution, and will be further discussed 
in Section| 6 T] where we show results. 



3.2 Deceleration in the adiabatic stage 

We use the lKompaneetsI ( Il96(]h approximation alongside the coor¬ 
dinate transformation identified by K92 to follow the adiabatic de¬ 
celeration of the shock front in a single power-law density profile. 


The assumptions and main equations of this prescription will also 
be used in our full numerical treatment for arbitrary density profiles 
(SectionID- We shall give here the analytic solutions for cu = 1 and 
3. These solutions will be used to validate our numerical treatment 
(SectionID- They also give an indication of the shock behaviour in 
a broken power-law density profile, when it expands fully interior 
or fully exterior to the density break. 

The Kompaneets approximation involves setting the post- 
shocl 0 pressure, P', to be uniform throughout the shock volume 
and equal to (some fraction. A, of) the mean interior energy den¬ 
sity. For an arbitrary volume V, 

(4) 

wherein the Kompaneets approximation proper is to take A to be 
constant. The ratio of specific heats is taken to be 7 = 5/3 both 
internal and external to the shock. 

Two additional assumptions in the treatment are that the direc¬ 
tions of the local velocity vectors along the shock front are normal 
to the shock front, and that the magnitude of the velocity is deter¬ 
mined by taking the post-shock pressure to be equal to that of the 
ram pressure of the environment (pWg, where p is the density of the 
unshocked gas) at that point (K92); 


Vs{R,t) 


(72 - 1) \E 

MR) V{t) ■ 


(5) 


Following the coordinate transformation of K92, the ‘time’ is 
parametrized by y (which actually has a dimension of length) via 


Ay 


( 7 ^ - 1) \E 
MV{y) 


( 6 ) 


as well as the dimensionless parameter x = \2 — uj\y/ ( 2 a) = 


® The three-dimensional volume integrals (of an offset sphere) over a sin¬ 
gular density converge for the shallow power laws used here: 1/2 ^ iUin ^ ^ For thermodynamic variables, we use primes (') to indicate the post- 

3/2 for p (X shock values (the values behind the shock front). 
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y/yc- The parameter x is, therefore, equal to y scaled with respect 
to a critical value t/c, which is when the shock either reaches the 
origin (lo = 1) or ‘blows out’ to infinity {lo = 3). Therefore, x (like 
y) can be considered to represent the ‘time’ in this transformation. 
The constant A ~ 1 is given by the difference in pressure behind the 
shock front relative to the ave rage pressure i nternal to the remnant, 
and in a power-law profile is l lShaDirdll979h 


(17-4a;) /9 

1 _ (9 _ • 


(V) 


In an ambient density with a single power-law form of equa¬ 
tion 13, the K92 transformation gives a self-similar solution to the 
Kompaneets equation (see equation^ 10 and 11 of K92) for an ex¬ 
plosion at i? = a (see Fig.[T](: 


( 


R 

a 


) 


2a 


- 2 



a 

COS [aO) 


— x^ + 1 = 0, 


( 8 ) 


for the polar coordinates R and 6, where a = (2 — a;) /2. This 
can be identified as a circular solution for a given x in the two 
variables {R/af and a9. Analytic solutions for the volume, time 
and velocity in a; = 1 and 3 densities are presented in Appendix IB] 
The equations describing the shock front can alternatively be 
parametrized by tpo, the initial angle of a point on the shock with 
respect to the axis of symmetry. The subsequent equations of mo¬ 
tion for a given tpo describe the paths of flowlines in the shock in 
terms of the polar coordinates measured from the SMBH (K92): 


2 \l/( 2 a) 


i? = a (l -b 2xcos'ipQ + x^) 
g - ± arctan ( 

U — 7- 7 dlC/Ldll I -— 

\a\ \1X cos'ipo 


(9) 

( 10 ) 


This flowline-based treatment is a useful context for the numerical 
approach to the shock evolution presented in Section|4A] and these 
equations will be used to compare with the numerical results. 


3.3 Intermediate-asymptotic transition 

The Kompaneets solution for the velocity diverges for a; —>■ 0, 
given that the volume V{x) —0. In this limit, the energy den¬ 
sity and, therefore, also the velocity, tend to infinity. The solution 
is, however, not intended to describe the initial evolutionary stage 
of the remnant. In order for the numerical treatment to correctly 
follow the SNR evolution, we must account for the initial coasting 
stage. A full analytic joining of the intermediate-asymptotic solu¬ 
tions between the ejecta-dominated and adiabatic stages is com- 

plex, even for an uj _= 0 ambient medium (see, for example, 

iTruelove & McKe3ll999h . 

As a model for this intermediate behaviour, we employ an ef¬ 
fective density (mass) term to the solution that gives a transition 
between the expected solutions. The density of the medium is mod¬ 
ified to: 

Peft = P(i7) + ^. ( 11 ) 

where the additional effective term counters the divergent be¬ 
haviour of the velocity at small volumes. This has the desired prop¬ 
erty that when the volume is large p — p{R) in the standard Kom¬ 
paneets approximation, while at small volumes the second term 
dominates to provide the initial coasting phase of the remnant. 


With this effective mass term, the SNR leaves the ejecta-dominated 
phase around the point at which the mass swept up from the envi¬ 
ronment is comparable to the initial mass behind the shock. 


3.4 Transition to the radiative stage 


As the shock slows, the late evolution of a typical SNR is marked 
by an increase in radiative losses. Although we will not model 
this stage, we intend to check the time-scales over which SNRs 
will reach this stage in quiescent nuclei (if they survive sufficiently 
long). 

Typically, cooling functions show a marked increase in ther¬ 
mal ra diation once the gas temperature drops to ~ 10® K (for ex¬ 
ample, ISchure et al.ll2009l) . This occurs due to the formation of a 
sufficient number of electrons bound to ions to allow for effective 
line cooling. Once regions of gas behind the shock drop to this tem¬ 
perature, the deceleration of the SNR becomes more pronounced. 
By calculating the temperature behind the shock we can determine 
the time at which parts of the remnant begin to cool more effec¬ 
tively. 

It is possible to determine the temperature of the shocked gas 
via the ideal gas law. 


P' = 


kBp'T' 

rriup 


( 12 ) 


(where, again, we denote post-shock values with primes, fee is 
Boltzmann’s constant, rriu is the atomic mass unit and p is the 
mean molecular mass), as well as the jump conditions for the (post¬ 
shock) density and pressure. 


; 7+1 

P = P - 7 

7-1 

For 7 = 5/3, 


and P' = 


‘2pVa 
7 + 1' 


(13) 


, 2 (7 - 1) 2 3mu/i 2 

- {j+lfkB “ 16 fcB 


(14) 


and the post-shock temperature is found to be T' ~ 10® K for Vb ~ 
300 kms”*^. Therefore, if we monitor each point along the shock 
for the time at which the velocity drops below this value, we may 
estimate the time at which radiative processes become significant. 

As the cooling function is also dependent on p, for a given 
temperature the rate of cooling is also expected to be amplified in 
regions of post-shock material with higher density. However, by the 
time that SNRs are radiative, they have survived the expansion past 
the SMBH and entered into the more uniform density beyond the 
SOI, such that the ambient density is similar across all points of the 
shock. At this stage, the SNRs are reasonably symmetric around 
the SMBH and the velocity is similar across all of the shock front, 
so that most of the SNR reaches the radiative stage at similar times. 
If the SNR survives expansion past the SMBH, this late evolution 
is largely uninfluenced by the details of any early interactions near 

the SMBH0. 


4 EVOLUTION OF REMNANTS AROUND QUIESCENT 
BLACK HOLES: NUMERICAL TREATMENT 

Purely analytic solutions for the shock front evolution via the Kom¬ 
paneets equation are not feasible for many density configurations. 


® Note that there are two sign errors in the exponents of equation 11 in ® For a detailed cons ideration of the radiative transition in power-law me- 

K92. dia, see lPetruj i2005li . 
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Figure 2. The basic numerical scheme, as described in Section ITTI Shaded 
boxes show the basic types of objects in the numerical constmction. The 
choice of the number of fiowlines determines the resolution, and only a few 
are shown here schematically. The flowlines track local physical properties 
of the shock (velocity and ejecta mass fraction). Collectively, they define the 
location of the shock front, with global physical properties such as its en¬ 
ergy density, determined by the volume. Along with the environment (most 
importantly, the background mass density), the global shock properties de¬ 
termine the evolution of all the individual flowlines in the subsequent time 
step. 

Therefore, we developed a numerical method that solves for the 
evolution of a shock front using the physical assumptions of the 
Kompaneets approximation described in Sectionj^ 

The primary assumptions that must be encompassed by the 
method culminate in constraints on the velocity. Namely, the di¬ 
rection of the velocity of any point must be perpendicular to the 
shock front, and the magnitude of the velocity must be determined 
by the energy density behind the shock and local ambient density 
as prescribed in equation 0- 

4.1 General prescription 

The numerical treatment follows an approach by which the shock is 
described by the evolution of flowlines through the background gas. 
The flowlines are the paths followed by tracer ‘particles’ (points) 
distributed along the shock front, analogous to the analytic treat¬ 
ment with the tpo parameter of equations © and oni- 

Fig. [2 shows a schematic of the approach. The initial (spheri¬ 
cal) state of the shock is broken down into flowlines characterized 
by their angle t/jQ. During the evolution, the number of flowlines is 
dynamic. To keep a reasonable resolution of the shock front, new 
flowlines can be inserted, with mean properties of adjacent flow¬ 
lines, if the distance between two points on the shock is over a 
defined threshold. For our simulations, a threshold of the order of 
0.05 pc has proven sufficient to describe a smooth shock front evo¬ 
lution on Milky Way-like scales. Flowlines may also be deleted 
in regions where parts of the shock front are colliding. The back¬ 
ground gas prescribes the evolution of the shock, but the behaviour 
of the post-shock gas is not tracked, and thus the background gas 
can be treated as being independent of the shock. 

The kinematics of the shock front are determined by the ve¬ 
locity vectors at each flowline. To determine the magnitude of the 
velocity, we use the jump conditions across the shock, and for those 


we need the energy density within the shocked volume (which is as¬ 
sumed to be a constant within this volume) as well as the local mass 
density in the environment (see equations Island [TTPI This stage of 
evolution is adiabatic, and so given an initial explosion energy we 
therefore calculate the energy density using the instantaneous vol¬ 
ume enveloped by the shock front. 

In an axisymmetric arrangement of gas density and explosion 
point, calculation of the volume is simplified by the geometrical 
symmetry; it is determined by a solid of rotation of the area of a 
two-dimensional slice about this axis. Any arbitrary ordered set of 
points {xi,yi) can specify the location of the shock front. Given 
these two-dimensional coordinates, the volume, by the second the¬ 
orem of Pappus, is equal to the product of the area of the non¬ 
intersecting polygon defined by these coordinates and the distance 
travelled by its centr oid under rotation about the symmetry axis 
jKem & Blanj[l94^ . Using the fact that the components of the 
centroid, C = (Cx, Cy), of a polygon are given by 

^ n — 1 

^ + Ci-ri) {xiVi+i - Xi+iyi) (15) 

i=l 

for ^ € {a;, y}, the volume of the SNR can be determined from 

n — 1 

V" = - ^ (t/i + yi+i) {xiyi+i - Xi+iyi). (16) 

i=l 

With the magnitude of the velocity known, each point on the 
shock evolves by determining the unit vector for the velocity that 
is perpendicular to its neighbouring points. The position is then 
linearly translated over a small time step using the velocity vector. 

We assign an ejecta mass element to each flowline. Fig. 
shows a schematic of this implementation. Due to rotational sym¬ 
metry, each point on the shock at f —^ 0 represents a ring segment 
of the SNR in three dimensions. We assign a thickness to each of 
these ring segments based on the spacing between the flowlines in 
the initial spherical state. The fraction of ejecta mass represented by 
the flowline is then the ratio of the area of this zone of the sphere 
to the total surface area of the sphere. 

Due to asymmetry in the background density and the presence 
of a strong density contrast near the origin, segments of the shock 
may collide with one another. This shock front self-interaction can 
lead to collisions in which kinetic energy is converted into inter¬ 
nal energy. Since it cannot be easily radiated away, we expect a 
transient acceleration outward of the heated gas, after which the 
fluid will return to the dynamics imposed by the global expansion. 
The numerical treatment of these self-interactions is outlined in Ap¬ 
pendix This treatment results in the deletion of some flowlines, 
accounting for the modification of the flow in this region. 

4.2 Comparison with analytic solutions for single power-law 
profiles 

Fig. a shows the morphology of an SNR, running into a circum- 
SMBH environment with density power-law gradient a; = 1 (left- 
hand panel) and a; = 3 (right-hand panel). There, we compare the 
analytic prescription described in Section lT^ and Appendix IB] with 
our numerical method. The numerical solutions are found to match 

^ For simplicity, the ratio A Ri 1 of the post-shock pressure to mean interior 
energy density is set to be exactly unity in equation 0. As the shock veloc¬ 
ity is proportional to y/X, the effect of this is small compared to other limi¬ 
tations inherent in the Kompaneets approximation discussed in Section l43] 


























Fate of SNRs near SMBHs 1 



Rx j Cl Rx / ex 


Figure 4. Numerical results (blue) compared with analytic solutions (red, dashed) for the locations of the shock front in density profiles with ui = I (left-hand 
panel) and a; = 3 (right-hand panel). The initial (spherical) state for the numerical solutions is shown in green. Units are given as ratios of distance {Rx, Ry) 
to explosion distance a from the density singularity at (0, 0) (the SMBH in our model). The results can be written in parametric form in terms of x, which 
increases with time t up to a critical value of a; = 1; see Section lX^ as well as the expressions for t{x) in Appendix I bI The solutions for a; = 1 are found up 
to x = 1, while for oj = 3 they are given up to x = 0.8 due to the divergence of solutions as x —t- 1 in this latter case. The trailing part of the shock (the 
part directly towards the SMBH) for tu = 1 solution reaches i?x/u = 0 at x = 1, while leading point (directly away from the SMBH) reaches i?x/u = 4. 
The trailing part of the ui = 3 solution asymptotically approaches a distance of R^/a = 1/4 as x —> 1, while in the same limit the leading paif of the shock 
follows Rx/a, —> oo. 



Figure 3. A schematic of the initial, spherical state of the SNR, with ini¬ 
tial positions for the flowlines (filled circles) around the point of explosion 
(open circle labelled ‘SN’). Flowlines in the positive-y portion of the x-y 
plane define a three-dimensional shock front by rotation about the axis of 
symmetry. Rotating each flowline about this axis produces a ring (shown 
for the fiowline as a thick line). The midpoints between the fiowline 
and its neighbours define the limits of the zone of the sphere assigned to 
that fiowline (thin lines). For a sphere of radius R ', the area of this zone is 
proportional to the height of the zone, h, since its surface area is 2nR'h. 
The fraction of total ejecta mass assigned to the fiowline is then the ratio of 
this area to the total area of the sphere. 

the analytic form very well. There is very slight deviation between 
the two methods, more noticeably in the w = 3 case, which is 
due to the fact that the numerical method requires a small spherical 
initial step. The analytic solution is closer to a sphere at small times 
in the uo = 1 solution so there is almost no discernible discrepancy. 

Fig.l^shows a comparison between our numerical (solid lines) 
and the analytic (dashed lines, arbitrary scaling) results in a broken 


power-law medium with cuin = 1 and oJout = 3. This figure shows 
the distance and velocity evolution of selected sample points on 
the shock front. We follow the portion of the shock that propagates 
towards the SMBH (blue lines), away from the SMBH (red lines) 
and at an initial angle of t/o = n/2 (green line). The numerical 
radius and velocity are seen to transition from the coasting (radius 
R' oc t, velocity v = const.) phase to forms similar to those seen 
in the pure Kompaneets solutions. 

As expected, the evolution of the trailing part of the shock, 
as it gets closer to the SMBH, approaches the analytic solution for 
a pure lj = 1 medium. Likewise, the leading part of the shock 
asymptotes to the pure a; = 3 analytic solution, as it expands 
away from the SMBH. The green line shows how the evolution 
of a fiowline that emerges at 90° from the 0 — 0 axis has, in¬ 
stead, an intermediate behaviour, which is influenced by the over¬ 
all broken power-law density. The figure also shows for reference 
the S edov-Taylor solution (black dashed line, lTavloilfl95(il : I Sedov! 
Il959l) for an explosion in a uniform ambient medium {R' oc 
and V oc 


4.3 Caveats and limitations of the model 

For more complex background density configurations, such as one 
with many large density contrasts that trigger self-interactions and 
turbulence, one may consider a treatment of self-interacting shocks 
that is more in-depth and sophisticated than that presented in Ap¬ 
pendix 13 The increase in velocity of any small self-intersecting 
region is expected to be a brief transient phenomenon; therefore, 
we do not presently apply any boost in velocity when merging 
flowlines, instead only accounting for the net direction of the flow 
that results from two colliding parts of the shock. The reason is 
that, in all our simulations, the portion of the shock front which 
undertakes self-interaction is limited, and therefore the treatment 
of these regions have a small effect on the overall volume evolu¬ 
tion. Obviously, if one considers a more complex geometry where 
self-interaction dominates the evolving volume, full hydrodynami- 
cal simulations are the only reliable tool of investigation. 
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Figure 5. Example of radius and velocity evolution for a remnant in a broken-power-law medium (solid lines). The explosion occurs at 3 pc (near the density 
break between cjin = 1 and cjout = 3) around a 5 x 10 ®Mq SMBH using the scaling described in SectionfJ] Solid lines are plotted from snapshots of the 
evolution of the remnant in the numerical treatment, where blue (lowermost) curves are for the trailing flowline (towards the SMBH) and red (uppermost) 
curves are for the leading one (away from the SMBH). The green curve shows the behaviour of a flowline in the numerical treatment that emerges at 90° from 
the 9 = 0 axis (-00 = ’n:/2). Some analytic results (with arbitrary scaling) are given as dashed lines for comparison. Black is the form of the Sedov-Taylor 
(uniform medium, = 0) solutions. The red dashed curve shows the solution for the point on the shock travelling directly away from the SMBH for a shock 
in a purely cj = 3 medium. The blue dashed curve shows the solution for the point on the shock travelling directly towards the SMBH in a purely u) = 1 
medium. 


The Kompaneets approximation itself has some drawbacks, 
in that it generally predicts too large a v elocity, and therefore 
size, for the shock on ce it accelerates i IKoo & McKed Il99d: 
iMatzner & McKe3ll999h . In the context of the present problem, 
this is more pronounced in the outer density region with a steeper, 
R~^, gradient. If much of the shock is in the R~^ region, the over¬ 
estimation of velocities and sizes will therefore be greater. 


5 GALACTIC NUCLEI MODEL 

In our quiescent SOIs, with no appreciable inflow of gaseous mate¬ 
rial from further out, the gas density distribution is that of an RIAR 
The distribution of early-type stars (the only population of interest 
here) is dictated only by local and current conditions, not bearing 
imprints of the long term history of the assembly of the nucleus. 
These facts will allow us to rescale features of our Galactic Centre 
(observationally constrained because of its proximity) to quiescent 
nuclei with different SMBH masses. 


of |Peeblesl ( ll972h . 


Rsoi = 


GM. 


(17) 


for a black hole of mass M,, where a is the velocity disper¬ 
sion of stars about the SMBH. We use this parameter not only 
to define the outer edge for the range of explosion distances 
considered, but also to rescale Milky Way properties to galac¬ 
tic nuclei with different Mm. To obtain an expression for the 
SOI which depends only on M,, we use the well-known (‘M,- 
a’) relation between black hole mass and v elocity dispersion 
dperrarese & Merrit3l200d : Gebhardt et al.ll200d) . Using the obser¬ 
vationally determined [Oebh^dFet^ ([20^ resultlfl M. = 1.2 X 


10® (ct/( 200 kms M©, we obtain 


^^soi 


M. 


4.3 X 106 Mq 


7/15 


pc, 


(18) 


where here, and hereafter, we rescale equations for the Galactic 
Centre black hole mass. For what follows, a useful parameter to 


5.1 Characteristic radii 

We consider SN explosions within the SOI of an SMBH. Their fate 
can be influenced by both the SMBH gravity and its gaseous envi¬ 
ronment. Correspondingly, there are characteristic radii in the nu¬ 
cleus associated with these properties. The first is that of the SOI: 
the range out to which the gravity of the SMBH dominates over that 
of the gravitational potential of the bulge. Following the definition 


® Recent studies imply that the AF,-cr relation is st eeper than this, 
and ( 7 may have an exponent closer to 5 (for example, iMorabito & Daa 
l2012h . Although there is still some ambiguity in the value of this ex¬ 
ponent, we tested the effect of a very steep relationship M, = 1.2 X 
10 ® (cr/(200 ® Mq motivated by IMorabito & d2012h . 

Even with this large exponent, we find that our main results, the time-scales 
in Section l631 are generally only increased by a factor of 2 (while the scal¬ 
ing of radii by the spheres of influence also increases by at most a factor of 
2). As the overall consequence is small, we do not present additional results 
for a steeper M.-cr relation in this work. 
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rescale the Milky Way properties is the ratio ( oc between 

the Rsoi of a generic M, and that of Sgr A*. 

The closer a supernova explodes to the SMBH, the stronger 
the tidal forces, which may become high enough to disturb and 
eventually disrupt the remnant in a dynamical time. This happens 
when the velocity of the shock front becomes comparable to the 
Keplerian velocity, wk, associated to the SMBH gravity field. We 
do not model distortions due to tidal effects but we account for 
the tidal disruption of the remnant when we quantify its ‘lifetime’ 
(see Section l631 l. To this end we test for whether Us < vk to detect 
parts of the shock that have decelerated enough to be sheared by the 
SMBH. We therefore introduce another characteristic radius—the 
innermost radius for the existence of SNRs, Rsh, which is limited 
by SMBH shearing. This minimal shearing radius is the point at 
which Vk is comparable to the initial SNR ejecta velocity, Uinit: 


the SOI. Since the total mass of stars in the SOI scales with M,, 
then so too will the number of massive stars and, therefore, the 
accreted mass: M oc M,. As above. Merit oc M,, as well, and 
therefore the ratio M/Mcrit is constant over M,. In other words, 
for our physically motivated picture of quiescent nuclei, the SMBH 
is accreting at the same fraction of Eddington as Sgr A* (M ~ 
lO-®MEdd)0 

Given this accretion rate of M/Mcrit ~ 10“®, we can esti¬ 
mate the density at Rq for nuclei with a different M,. To do so we 
use the continuity equation for the flow, 

M « 4nRl p{Ro) vk{Ro), ( 21 ) 

where the scale height H fa R and the radial velocity ur ~ uk- 
The reference density for our general galactic nuclei models 
is therefore 


Rah = 


GM. 


= 1.9 X 10” 


( ^ 

/ Pinit \ 

V4.3 X lOSMoy 

VlO^ kms“M 


= 900 ( ^ ^ Rfr, 

VlO^kms-V 


(19) 


where R^ is the gravitational radius of the SMBH. 

Supemovae that occur within Rsh are completely sheared, as 
the velocity of the ejecta in all directions is less than the Keple¬ 
rian velocity around the SMBH. Note, however, that SNRs can be 
sheared also at larger radii as the ejecta slows down, and may reach 
the local Keplerian velocity at a radius larger than Rsh - Comparing 
the shock velocity to the Keplerian velocity is effectively equivalent 
to comparing the ram pressure with the ambient baryonic pressure, 
Pgas, since uk ~ Cs oc ^/P^IsJp (where Cs is the sound speed 
in the external medium). Additionally, we note that for the same 
reason we can ignore the shearing of remnants during the initial 
explosion for a > Rsh, we can also neglect the (Keplerian) orbital 
motion of the progenitor stars. 

Finally, for the gas models, the reference radius for the density, 
Ro, and the location of the break in the gas density power law Rh 
(as explained in Fig. [TJ are scaled in our model by the SOI, such 
that Ro = CRo,mw and Rh = C^b.MW- 


5.2 Gas models 


As mentioned previously, we expect that quiescent SMBHs are sur¬ 
rounded by RIAFs, similar to that which is suggested in the Galac¬ 
tic Centre. This implies that all flows have similar density gradients, 
which is set by the physical processes which characterize this accre¬ 
tion regime. Additionally, their accretion rate must be modest and, 
in particular, lower than the critical value for advection-dominated 
accretion of Merit = a^Msdd, or 


Merit = 9 X 10 


M. 


4.3 X 106 Mq 


Mq yr 


( 20 ) 


where Msdd is the Eddington rate, a = 0.3 (iNaravan & Yill 199^1 
and we assume a 10% radiation efficiency for MEdd- At a reference 
distance of 7?o,mw ~ 0.04 pc, the accretion rate is estimated from 
observations and simulations to be around M « 10“® M© yr~^ 
jCuadra et al.ll200d : lYuanll2007ll . Therefore, Sgr A* is accreting at 
~ 10 ~® of its critical rate. 

We extend properties of the Sgr A* accretion flow to other 
quiescent nuclei as follows. The primary material for accretion in 
quiescent nuclei originates from the winds from massive stars in 


"° = "(^°)"130( 4.3x^06Me ) ^^2) 

An example of our density model for M, = IO^Mq and three 
different inner gradients is given in Fig.|^ where the effect on the 
density profiles of fixing the scaling reference point at Ro is evi¬ 
dent. 


5.3 Massive star distributions 


Given a physical number density nt,{R) of stars at a distance R 
from the centre of mass, the projection on to the celestial sphere 
gives, as a function of the projected radius i?pr, a surface density 
of stars E*(i?pr). Observationally, the latter quantity is typically 
given. Assuming a spherically symmetric spatial distribution, it is 
possible to reverse the projection to infer the spherical number den¬ 
sity. 


n*(7?) 


^ dS« {Rpf dRpr 
^ Jr di?p. - R^ ’ 


(23) 


provided that the physical number density n* {R) falls off at large 
J? at a rate greater than R~^.For power-law distributions, this gives 
a correspondence of the observed radial dependence, RpJ^, to the 
physical dependence, R via the relationship 7 ~ T -F 1. 

Eor the Milky Way, there is evidence for two different power- 
law distributions in the old and young stellar populations of the 
Galactic Centre. A particular curiosity is an apparent deplet ion of 
late-type (K, M) giants in the inner 0.5 pc i IDo et alj|20d^ . This 
leads to a much shallower (possibly inverted) inner power-law for 
the late-type distribution compared to that of the early-type (O, 
B) stars. Recent analyses estimate the radial dependence of the 
early-type stars in the Milky Way nuclear star cluster to be ap- 

S iroximately Rttr ins i de the power-la w break and Rpf^ outside 
Buchholz et aHl2009l: lDo et al.l2013al) . corresponding to values of 
7 of 2 and 4.5, respectively. 

In the nuclei of different galaxy types, variations in the stel¬ 
lar distributions for longer-lived stars are possible due to differ¬ 
ing nuclear assembly histories. However, in our picture of a self¬ 
regulating SOI, the young star distributions are taken to be the same 
across the range of M,, where the most recent star formation in this 


® It is possible for SMBHs to be accreting at different fractions 
M /Merit < 1 and still be termed ‘quiescent’ in the conventional sense. 
However, the accretion rate given here is the most physically motivated 
value based on scaling of quantities by M,, and deviations from this value 
are beyond the scope of this work. 
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Figure 6 . Example of a gas density model used for a galaxy with a IO^Mq 
SMBH, showing the number density, n, as a function of radius R from 
the SMBH. The red (shallowest) line corresponds to an inner gradient with 
= 1/2, the green to ojin = 1 and the blue (steepest) to cjin = 3/2. 
All models have a gradient outside the break of cjout = 3. The density is 
scaled using a reference point Rq, seen as the point of convergence of all 
the inner density gradients (in this case, Rq = 0.06 pc). A break in the 
density distribution is located at a constant R = for all choices of the 
density gradient (in this case, R\^ = 0.6 pc). The left-hand and right-hand 
limits of the horizontal axis are determined by the shearing radius (equation 
[19) and SOI (equationflT^. respectively. 

region is indifferent to the history of the nucleus. Therefore, for our 
galactic nuclei model, we use the same values for 7 as those given 
above for the early-type stars around Sgr A*. As before, we scale 
the break in the stellar number density by the SOI of the SMBH, 
which defines the transition radius between the two values of 7. 


6 RESULTS 

We proceed to describe our main results, based on the method out¬ 
lined in the previous sections. In Section Ibdl we examine the ef¬ 
fect of black hole mass and gas density profile on the deceleration 
length using the prescriptions of Section [JTI Then, using the nu¬ 
merical method of Section |4] the overall SNR morphology is pre¬ 
sented in Section lb^ Finally, in Section l63] we investigate the X- 
ray emitting lifetimes based on shearing of the SNR ejecta by the 
SMBH. We then use this last result to predict the mean SNR life¬ 
times expected within the SMBH spheres of influence for different 
M,. 

6.1 Deceleration lengths 

The SNR begins to appreciably decelerate once the swept up mass 
becomes comparable to the ejecta mass. This end of the ejecta- 
dominated stage can be characterized by a deceleration length from 
the explosion point, L, determined by the density integrals of Sec¬ 
tion [TT] which varies with direction. This deceleration length de¬ 
pends on the gas density and on the radial density profile. These de¬ 


pend respectively on the SMBH accretion rate and accretion mode 
(for example, CDAF versus ADAF). 

Considering the values of L in various environments, we 
obtain an indication of the length- and time-scales over which 
SNRs will end their ejecta-dominated stage and start decelerat¬ 
ing. As a reference for the crossing time-scale of a nucleus, re¬ 
call that an SNR that does not appreciably decelerate from its ini¬ 
tial ~ 10"^ kms“^ would reach a radius of 1 pc in approximately 
100 yr. An investigation of L in different directions indicates which 
SNRs will decelerate within this time-scale. It also provides a test 
for the level of asymmetry of the SNR during this stage of evolu¬ 
tion. We will later proceed to model SNRs through the decelerating 
stage. 

Figs □ and [8] depict two curves describing the deceleration 
length, L, approximately towards and away from the SMBH. The 
angle tpo approximately towards the SMBH is taken to be 10“® tt, 
such that the integrated path through the density runs very close to 
the SMBH, but does not pass though the singularity at the origin. 
The difference between these two curves provides a measure of the 
asymmetry of the remnant at the end of the ejecta-dominated stage. 
For comparison, in Fig. |7] a third curve (solid line) is shown that 
describes an average deceleration length derived by an integral over 
a sphere. 

We first consider a model of the background gas in the Milky 
Way. Fig.|7]shows L for three different density values: the observa- 
tionally motivated ‘canonical’ density po = pc (in number density, 
ric ~ 130 cm“®, left-hand panel), and 3 (central panel) and 30 
(right panel) times that value. We consider values other than the 
canonical density, as, even in quiescent nuclei such as the Galactic 
Centre, there is the possibility for variation in the overall density 
of the accretion flow. For example, denser accretion flows can re¬ 
sult from sudden accretion episodes from tidally disrupted stars or 
cloudsl3 or they can be associated with more intense star forma¬ 
tion activity in the nucleus. Scaling the density also shows the effect 
of under- or misestimating the gas density from the X-ray emission. 

As expected with increasing density, there is an overall trend 
towards lower values of L/a. There is also a trend towards more 
symmetric remnants with increasing density, since in general the 
ratio L/a is reduced for higher densities. 

The investigation of different density profiles (see Fig. [Hi leads 
us to conclude that C DAF/ADIOS mod el, preferred by Galac¬ 
tic Centre observations dWang et al.ll2013h . gives, quite generally, 
shorter deceleration lengths (red lines in Figures|6]and|7). The flat¬ 
ter CDAF/ADIOS profile (smaller ojin) is denser in most of the SOI 
of the black hole, therefore reducing L/a. 

For the canonical value of density in the Milky Way (left panel 
of Fig. |7), the deceleration lengths are L > 1 pc. Considering the 
CDAF model, we remark that, for the canonical density, the ma¬ 
jority of the SNRs would decelerate beyond the SMBH location. 
Ejecta from a star such as S2 (marked with a blue dot-dashed line 
in Fig.jT} is expected to evolve more symmetrically than that from a 
star further out, in the stellar disc(s) (~ Rh). Already with a factor 
of few enhancement in density, SNRs in and beyond the stellar disc 
would decelerate appreciably before they reach Sgr A* (see central 
panel). 

Fig. in shows deceleration lengths for each of Uin € 
{1/2, 1, 3/2} for galactic nuclei with SMBH masses of 10^ M© 


A recent exam ple around Sgr A* is the object G2 (for example, 
iBurkert et al.ll2012|) ; though, if a cloud, its mass is too small to have a sig¬ 
nificant impact on the overall density. 
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Figure 7. Deceleration lengths, L, as a function of explosion distance from the SMBH, a, for a Milky Way model (M, = 4.3 X 10®M@) with density 
scalings of the canonical density 130 cm~® of 1, 3 and 30 times (left to right). Red lines show the 0 ;;^ = 1/2 model and green lines show the 0 ;;^ = 1 
model. The solid curves are the average deceleration lengths derived by integrating the mass over a sphere, as discussed in Section iTT] The dotted and dashed 
curves are integrals away from and towards the origin, respectively. The shearing radius (2 X 10“^ pc) and the SOI for Sgr A* are marked by and i?soi, 
respectively. The vertical dotted black line is the location of the break in density, outside of which tUout = 3. The diagonal dashed black line shows L = a, 
which represents the distance to the SMBH. The blue dot-dashed line shows the pericentre distance of the star S2. The curves for tUin = 1 are not plotted in 
the left-hand panel since most of the points lie outside of the SOI axis bounds. 


(red), 10® M 0 (green) and 10® Mq (blue). Hereafter, we scale the 
explosion distance by the location of the break in gas density, R^, 
as it is the value of a at which a change in behaviour of SNRs is 
expected. Accretion rates (and thus gas density normalization) and 
characteristic radii are rescaled as explained in Section and all 
increase with black hole mass. 

For M, < 10® M 0 , the density is low enough and 7?ah is 
small enough that for small values of a the SNR can expand over 
the SMBH before appreciably decelerating. The centre of such a 
shock front is close to being aligned with the centre of symmetry 
of the gas distribution, which leads to a more symmetric evolu¬ 
tion. This is seen in the fact that all L values converge at small a. 
However, by M, > 10® M 0 , the gas density is high enough and 
Rsh extends far enough from the SMBH that explosions near the 
SMBH cannot pass over the singularity before being sheared, and 
this more symmetric expansion regime beyond the SMBH at small 
a is no longer present. 

In general, with increasing overall density (around more mas¬ 
sive black holes), we find the same trend observed for the Milky 
Way with denser gaseous environments: the ratio of Lja decreases, 
as does the maximum possible asymmetry (differences between up¬ 
per and lower curves). On the other hand, increasing ujin also tends 
to create greater asymmetries in the SNRs, as can be appreciated 
by comparing the three panels of Fig.[^ For higher tJin, the ratio 
L/a is higher in the direction away from the SMBH as a result of 
the density being lower at any point further than i?o, but it is lower 
towards the SMBH as the SNR sweeps through a steeper density 
gradient near the origin. 

6.2 Morphological evolution 

We turn now, using the numerical treatment of Section [T] to the 
subsequent adiabatically decelerating evolution of the remnant. 
We consider explosions both inside and outside the density break, 
Rh- The explosion distances are chosen such that, across all M., 


the same a/Rh ratio is maintained for the examples inside i?b 
(o/i?b = 0.6) and outside i?b {a/Rh = 2.5); this will also be 
useful for comparison with the time-scale plots of Fig. m 

Fig. [ 9 ] depicts the morphology for two explosion distances in 
the Milky Way environment. The explosion inside the density break 
does not decelerate before reaching the SMBH (as expected from 
Fig.lll where L > a). Therefore it passes over the SMBH without 
any significant distortion. The SNR subsequently expands into the 
lower density region almost spherically, with the centre of the SNR 
being very near the SMBH. For the explosion outside the break in 
density, the trailing part of the shock decelerates before reaching 
the SMBH, while at the same time the parts expanding through the 
a; = 3 region wrap around the SMBH and eventually self-interact. 
SNRs such as this, which explode far enough from the SMBH that 
L < a, show significant asymmetries during their evolution. 

Fig-Hnishows examples of the variation in morphology of the 
remnant arising from differences in the black hole mass (and there¬ 
fore the gas density). The values of M, are the same as those used 
in Fig.[^for the deceleration lengths; the black hole mass increases 
from top to bottom. The data shown in Fig. [T^ are summarized in 
Tables 

It is evident that, in many cases, much of the mass of the rem¬ 
nant remains near the SMBH due to the focusing effect of the den¬ 
sity gradient on the flowlines. Therefore, unlike expansion away 
from the SMBH, where the mass behind the shock is more tenu¬ 
ous due to the rapid shock expansion, the ejecta material near the 
SMBH is expected to be more concentrated. 

In addition to the symmetries in L found in Section |6d] ex¬ 
plosions closer to a lower-mass SMBH (M, < 10® M©) are also 
found to be more symmetric during their adiabatic evolution com¬ 
pared to those further from the SMBH (compare panel ‘a’ of Fig. Ei 
to panel ‘b’). SNRs near lower M, expand over the SMBH largely 
unimpeded, and their centres are closely aligned with the centre 
of symmetry of the gas distribution when they enter the adiabatic 
phase. 
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U^in — 1/2 ^in — 1 




Win — 3/2 



Figure 8 . Deceleration lengths, L, as a function of the ratio of the explosion distance, a, to the break in gas density, R^. Each set of axes is to the same scale, 
and corresponds to a different gradient near the SMBH: cjjn = 1/2 (top left), = 1 (top right) and cjjn = 3/2 (bottom). As indicated in the top-left 
panel, red corresponds to the gaseous environment of an SMBH of 10^ Mq, green to 10® Mq and blue to 10^ Mq. Two curves are shown for each M* 
(each colour); the higher curve shows L for the direction away from the SMBH, and the lower curve shows L towards the SMBH. The dashed coloured lines 
correspond to i?sh value of Af*. The dotted black line shows the break in gas density at while the solid black line shows the extent of the SOI 

(which is the same multiple of i?b for all M*). Unlike in Fig. [7] we also show results for a standard ADAF model (cjin = 3/2; bottom panel), since these 
density profiles can produce L values that fall within Rsoi- 
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a/Rh 0.6 a/Rh ^ 2.5 




Rx [pc] Rx [pc] 

Figure 9. Morphological evolution for two explosion distances for the Milky Way, inside (left) and outside (right) Rj,. The explosion distance in the left-hand 
panel is a = 0.24 pc, and in the right-hand panel, a = 1.0 pc. The black hole is marked with a black point at the origin and the supernova occurs at the green 
star. The break in density at R\^ is shown as a thick blue circle. Snapshots of the shock front at different times are in solid black, and dashed lines show sample 
flowline paths. Flowlines flagged as having reached the shearing condition discussed in Section 1631 facMp < i>k) are shown as red points. Some example 
snapshot times are indicated next to corresponding dot-dashed magenta lines. In the left-hand panel, the initial snapshot shown is at 13 yr and the final at 
620 yr; in the right panel, the initial snapshot is at 13 yr and the final at 1100 yr. The spacing between snapshots is an equal multiple of the (variable) time 
step used. In the left-hand panel, the initial snapshot spacing is 16 yr, and the final spacing is 84 yr; in the right-hand panel, the initial snapshot spacing is 
16 yr, and the final spacing is 120 yr. In both cases, the remnants expand well beyond the window shown, until they reach the radiative stage as discussed in 

Section r6.3.1l 


panel 

M. (Mq) 

a (pc) 

finit (yr) 

Atinit (yr) 

Afflnal (yr) 

tfinal (y^) 

a 

1 X lO’’ 

0.28 

14 

23 

160 

1200 t 

b 

1 X 10’^ 

1.2 

28 

32 

1000 

3300 t 

c 

1 X 10® 

0.83 

19 

91 

760 

4300 

d 

1 X 10® 

3.5 

24 

58 

580 

5100 

e 

1 X 10® 

2.4 

16 

40 

210 

1000 

f 

1 X 10® 

10 

70 

170 

900 

4300 


Table 1. Parameters used for the plots of Fig. 1101 where M, is the mass of the SMBH and a is the distance of the explosion point from the SMBH. The 
snapshots that are plotted are at fixed multiples of the time step, although the time step is dynamic and increases with the size of the remnant. The first curve 
that is plotted is at a time iinit, and the difference between the first two snapshots is given in Ai. The difference between the final two snapshots is given 
in Aifinai, and the last curve that is plotted is at a time tfinal- Additional examples of snapshot times are shown on Fig. llOl with corresponding magenta 
dot-dashed curves. Note that the values of At are not the time step used by the numerical treatment, which is much smaller, f; simulation ran until 10“^ yr, 
stopping due to the radiative onset threshold discussed in Section 1631 but the last curve shown within these axis limits is at the stated time. 


Remnants around higher-mass SMBHs (M, > 10® M©) are 
more symmetric for the whole duration of their adiabatic lifetimes. 
The overall increase in density causes the scale of the remnant to 
be small (relative to the scale of the background gas distribution), 
and so, within the lifetime of the remnant, significant asymmetries 
have not yet developed. Indeed, by looking at the overall remnant 
size with respect to M,, there is a clear trend towards decreasing 
remnant size with increasing M, during this phase of evolution. 


6.3 Adiabatic SNR lifetimes 

During the adiabatic phase, the SNR can be observed as a hot, X- 
ray emitting object. The hard X-rays can penetrate the obscuring 
matter in galactic nuclei and allow us to detect the SNR. In partic¬ 
ular, in quiescent nuclei, the SMBH light may not prevent the SNR 


detection in X-rays. We will refer to this temporal window during 
which the SNR can be observed in X-rays as the ejecta ‘adiabatic 
lifetime’ or simply its ‘lifetime’, because once the ejecta become 
radiative it cools rapidly and its X-ray emission ceases. 

As previously mentioned, in a general environment, the adia¬ 
batic phase ends when the expansion has caused the internal tem¬ 
perature to decrease enough for radiative losses to become dynam¬ 
ically important. In the SOI, however, the tidal field can tear apart 
the ejecta well before the end of its radiative phase. In this case, the 
ejecta is dispersed and may not be immediately identified as such, 
at any wavelengths. 

So far, we have been evolving the shock expansion without 
directly considering the gravitational field of the SMBH, but only 
considering its indirect influence in shaping the gas density profiles. 
This is a very good approximation as long as the internal pressure 
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a/i?b ~ 0.6 a/Rh ~ 2.5 








Figure 10. Example remnant morphologies for three black hole masses in an a;in = 1/2 density for explosions at two a/i?b ratios (also indicated as dotted 
lines on Fig. ED. All axes are in units of parsecs. The top row (panels ‘a’ and ‘b’) shows M* of 10^ Mq, the middle (‘c’ and ‘d’) 10^ Mq and the bottom 
(‘e’ and ‘f’) 10^ Mq. Specific details of the parameters in each panel are given in Tabled All markers, line styles and colours are as for Fig.f^ Snapshots are 
plotted up until the end of the life of the remnant discussed in Section [63] fif found). 
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forces are significantly larger than the SMBH tidal forces. From 
momentum conservation, this means that we can ignore the SMB FI 
gravitational field whenever the shock front is faster than the local 
Keplerian velocity, vk- When, however, Va < vk, the dynamics of 
the shock front is dictated by gravity. As different adjacent parts 
move at different speeds ~ vk oc the SNR is sheared in a 

dynamical time. 

If the overall density is large enough, the swept up mass causes 
the shock front to decelerate before being sheared. As shown in 
Sections |6. II and [6^ most of the ejected mass is focused close to 
the SMBH (see Fig. [Toll, in regions of high Keplerian velocities, 
where this is more likely to occur. This effect is more prominent in 
the surrounding of more massive black holes. 

Practically, we consider that the SNR has ended its life when 
half of the original ejecta mass satisfies the condition Va < uk. 
A fraction of the ejecta mass is assigned to each flowline (as per 
Section lrTt . and we monitor the total proportion of sheared mass 
Msh to determine this time. The sheared portions of the remnants 
are marked with red points over the flowlines in Fig. m Fig.im 
describes the variation in lifetime for different explosion distances 
and M,. For models with M, < IO^Mq, no SNR within the SOI 
ends its lifetime by shearing within the 10^ yr shown on the plot. 

There is a clear trend in the behaviour of the remnant life¬ 
times as a increases. The lifetimes for SNRs exploding within Rah 
vanish, by the definition of the shearing radius. Just outside Rah 
the lifetimes are short but they increase rapidly with a. This is be¬ 
cause, at small a, the SNR expands at Vinit ~ 10^ kms“^ for the 
time that it is near the SMBH. Therefore, the amount of sheared 
ejecta is directly proportional to the fraction of the surface area of 
the SNR that enters the sphere of radius Rah around the SMBH. 
This quickly decreases as a increases, causing the rapid increase in 
lifetime. 

For M, < 10® M©, there is a sudden jump to very long life¬ 
times at distances slightly larger than Rah (in fact, within a range 
of small a, the Msh ^ Mej /2 condition is never met before the ra¬ 
diative stage sets in). Explosions at small a expand over the SMBH 
without any significant disruption due to their high initial veloc¬ 
ity. They then almost entirely travel ‘downhill’ in density, and so 
^ 50% of the remnant is never sufficiently decelerated by travel¬ 
ling into a sufficiently high densitvF^ 

Further increasing a, the lifetime then drops significantly. This 
can be attributed to the aforementioned focusing of ejecta towards 
the SMBH (as seen in Fig.llOt. Since the shock decelerates signif¬ 
icantly in these same regions, a drop in the lifetime is seen, partic¬ 
ularly near the break in gas density, Rh, in Fig. m For lower-mass 
SMBHs like Sgr A* (M. < lO'^ Mg), although 50% of the SNR 
is not found to be sheared, parts of SNRs at this region of a may be 
decelerated enough before reaching the SMBH that at least a small 
fraction (< 20%) of the remnant is sheared. At large a values, much 
of the SNR spreads out to distances further from the SMBH before 
appreciably decelerating. Here, a combination of lower vk and less 
deceleration of the shock conspire to lengthen the time taken to 
reach the shearing condition. 

Inspection of the left-middle panel of Fig. shows an inter¬ 
mediate behaviour around M, « 10® Mg . In this case, for explo¬ 
sions near the SMBH, the SNR can first pass beyond the SMBH 
and then become significantly sheared. However, once the SMBH 
mass is high enough (M. > 5 x 10® Mg), the higher densities and 


The amount of deceleration in such cases may be underpredicted by the 
Kompaneets approximation, as discussed in Section|4] 


larger vk mean that it is impossible to have large lifetimes for small 
a, as remnants are always significantly sheared before they expand 
far from the SMBH. As seen in Sections 16.11 and 16.21 they there¬ 
fore never enter the more symmetric, ‘downhill’ expansion regime 
seen for lower values of M,. Therefore, around these higher mass 
SMBHs, the lifetime is generally ~ 10® yr or less in the entire 
inner density region. 

We note that the most likely distances for core-collapse super¬ 
novae to occur are at or within Rh, as the density of massive stars 
falls much more steeply beyond this radius (see Section [6.3.2b . The 
SNR lifetimes are shortest in this region (see Fig. [ID- For more 
massive SMBHs the lifetime is also reduced significantly over the 
whole inner region. 

Before concluding this section we note that the tidal shearing 
of an SNR might have a significant observational signature. One 
can expect that the sheared ejecta will be largely accreted by the 
SMBH. This will temporarily enhance the accretion rate, leading to 
a period of about 100 years in which a few solar masses are accreted 
by the black hole. This is a very large accretion rate that may lead 
to a flaring of the SMBH, reaching luminosities comparable to the 
Eddington luminosity over this period. 


6.3.1 The radiative transition 

As we are interested in the hot. X-ray stage of evolution, we test 
for the onset of increased radiative losses in the remnant. This 
was estimated by comparing the shock velocity to the threshold 
of 300 kms“® outlined in Section [3^ As with the shearing con¬ 
dition, we consider the adiabatic lifetime to have ended when more 
than half of the remnant is below this velocity. 

For SNRs that are not destroyed by shearing, we find that this 
criterion is not satisfied before the SNR grows larger than the SOI. 
However, the density gradient will not continue to have a i?~® form 
indefinitely. Beyond the SOI, the SMBH no longer has any sub¬ 
stantial influence on the environment and the density is expected to 
level off. 

Therefore, to determine the onset of the radiative stage, we 
find the time when more than half of the shock reaches 300 km 
by extending the density uniformly beyond the SOI. If the SNR sur¬ 
vives to these large radii, we find that they eventually expand nearly 
spherically, and that the shape in this late stage is largely indiffer¬ 
ent to the processes that occurred near the SMBH. With a uniform 
density and approximately spherical evolution, the late-time kine¬ 
matics closely follow the Sedov-Taylor solution. For densities that 
flatten outside the SOI to about n ~ 1 cm~®, we find that the SNRs 
transition to the radiative stage at ~ 10"^ yr, independently of M,. 
All adiabatic lifetimes therefore end at this age if the SNR has not 
already been destroyed by shearing from the SMBH. 


6.3.2 Mean lifetimes 


Deriving properties over the whole SOI is useful in the cases where, 
for more distant galaxies, individual SNRs may not be observation- 
ally resolved. Young SNRs will contribute to the total X-ray emis¬ 
sion from these regions, which can be more easily observed. 

To summarize the effect of the SMBH environment on SNRs, 
we calculate the lifetime for core-collapse SNRs averaged over the 
entire SOI: 


(fad) 


ncc{R) R^ dR 
nac{R) R^ dR 


(24) 
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Figure 11. Time taken to shear a total ejecta mass M^h. ^ ^ej /2 for a range of M* and cjin, for explosion distances a scaled to the break in gas density, i?b • 
Decreasing tad curves correspond to increasing M%, and the black hole masses shown are as follows; magenta: 5 x 10^ Mq, yellow: 1 x 10® Mq, cyan: 
5 X 10® Mq, red: 1 x 10® Mq, green: 5 x 10® Mq, blue: 1 x 10^® Mq. The curve for each M» represents data from simulations covering between 30 
and 40 values of a. Dashed, vertical coloured lines show the con'esponding Rsh - The dashed black line shows the location of /?b- The dotted black lines show 
the two regions of explosion distances investigated in Fig. Qo] The right-hand limit of the horizontal axis is set at the SOI. By 10“^ yr, radiative losses become 
significant in all cases; if the shearing condition i;sh < has not yet been met by this time, the adiabatic stage ends due to the radiative transition discussed 
in Section f6.3.1l Note that the Galactic Centre SMBH is not shown on this plot, as we do not find Mgh ^ Mej/2 at any value of a before the radiative stage 
sets in. 


where 


ncc{R) oc < 



R Rb 
R > -Rb 


(25) 


is the volume number density of stars with mass over 8Mq (see 
Section|53J. 

The average lifetime, (fad) is shown in Fig.[T3 as a function of 
Mm. Examining Fig.[TT]and equation l l25b . we see that much of the 
reduction in (fad) is determined by the value of fad for explosions 
near and inside the break in gas density (a < Rb). The weighted 
contribution of fad to the mean lifetime is higher inside the break 
due to the higher densittF^ of massive stars; additionally, much 
of the reduction in fad occurs near the break for lower values of 
Mm, while fad is low throughout the inner region for higher M,. 
For Mm > IO^Mq, the mean adiabatic lifetime of the SNR gets 
increasingly shorter, well below the canonical value of ~ 10^ yr. 
By Mm > 1O®M0, the lifetime of most SNRs in the SOI is ended 
by disruption by the SMBH while the SNR is in the adiabatic stage. 
Shallower inner gas density profiles (green and red lines) amplify 
these trends. 


7 DISCUSSION AND CONCLUSIONS 

In this paper, we presented a novel numerical method based on the 
Kompaneets approximation for calculating the evolution of a shock 
in an arbitrary axisymmetric configuration of density. Our approach 
has the benefit of being more flexible than analytic solutions (which 



Figure 12. Estimated mean adiabatic lifetimes as a function of black hole 
mass in the present models, using equation (24) with the results for t{R) 
in Fig. llll and the stellar number density distributions n{R) in Section [531 
The red (lowermost) line is for oJin = 1/2, green for = land blue (up¬ 
permost) for a;in = 3/2. For lower M* (< 10^ Mq), remnants throughout 
the SOI are not destroyed by the shearing condition by the time radiative 
losses start to become significant. The adiabatic lifetime therefore ends at 
10 ^ yr due to the radiative stage, which is indicated by the dotted black 
line. 


Note that dN{M > 8) = ncc(-R)-R^ is constant for R ^ R\^. 


































only exist for some simple density stratifications) while being much 
faster than full hydrodynamical simulations. 

We apply this numerical method to trace the evolution of 
SNRs in quiescent galactic nuclei, with properties similar to those 
of our Galactic Centre. We describe these nuclei as self-regulating 
and steady-state systems, where gas inflow from outside this region 
has been limited and unimportant during at least the last 100 Myr. 
In this scenario, most of the star formation occurs in situ, recycling 
the gas ejected in supernovae, and winds from massive stars feed 
the accretion flow, which forms most of the interstellar gas. We 
predict the morphological evolution of SNRs in these nuclei, and 
relevant time-scales such as their X-ray lifetimes. 

We find that the supernova remnants that explode very near 
low-mass SMBHs (< 10® Mq, such as Sgr A*) will pass over the 
SMBH before appreciably decelerating and will continue expand¬ 
ing almost spherically. Although there can be prominent distortion 
during the early evolution of SNRs that explode further away from 
these SMBHs, SNRs at all explosion distances appear reasonably 
spherical by the onset of the radiative regime, much like SNRs in 
a typical interstellar medium. Notably, in the Galactic Centre, this 
implies that an SNR should be observable in X-ray for ~ lO'* yr. 
The pre sence of a suspecte d SNR enveloping Sgr A*, known as Sgr 
A East jMaeda et alj|2002h . fits with our prediction of SNRs being 
able to survive their adiabatic expansion through the Galactic Cen¬ 
tre region. 

If the SMBH mass is large enough (> 10® Mq), we instead 
expect a wide range of SNR morphologies, depending on the ex¬ 
plosion distance from the SMBH. There, SNRs typically end their 
life due to tidal shearing and disruption. The observable lifetime is 
therefore suppressed (10^ ~ 10® yr) with respect to SNRs evolv¬ 
ing around lower mass SMBHs. The reductions in SNR lifetime de¬ 
pend on the inner gradients of gas density as predicted by accretion 
theory. Conversely, therefore, observations of SNRs can be used 
to infer and constrain properties of their environment. For exam¬ 
ple, variations in the density gradients can produce different global 
quantities such as the mean lifetime of SNRs (Fig. m or their over¬ 
all sizes in their early evolution (Figs|3and[8}. 

The disruption of the SNR, or fraction of it, by the cen¬ 
tral SMBH, that takes place when this material slows down be¬ 
low the Keplerian velocity would lead to a period of enhanced 
accretion on to the central black hole. If a significant fraction of 
the sheared material is trapped by the SMBH we expect accre¬ 
tion of a few solar masses on to the SMBH over a period of 
~ 100 yr, yielding (assuming efficiency of 0.1) a luminosity of 
the order of 0.7 x 10“*^ ergs~^. This enhanced accretion would 
be at the sub-Eddington rate for the higher mass black holes (above 
0.5 X 10® Mq) but still very significant and at a level comparable to 
a powerful AGN. This may lead to a period of flaring of the other¬ 
wise quiescent black hole. Such events would happen even around 
the Galactic Center and other small-mass SMBHs. While SNRs 
are not completely disrupted around such black holes we still ex¬ 
pect events in which up to 20% of the SNR material is accreted 
on to the central black hole over a period of about a hundred years 
leading to significant flaring. 

Beyond the Milky Way, an excellent example of an SNR that is 
resolved in a galactic nucleus is S Andromedae (SN 1885A), which 
has an angular diameter of about 0.7 arcse c and has a morphology 
resolvable by the Hubbl e Space Telescope l lChevalier & Plaid 19881 : 
iFesen et al.l l999l , l2007h . The SNR is only 60 parsecs from the cen¬ 
tre of the bulge of the Andromeda galaxy, though not quite within 
the SMBH SOI. Although SNRs such as S Andromedae are resolv¬ 
able in other galactic nuclei with the current generation of instru- 
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ments, individual SNRs may not be distinguished in more distant 
galaxies. 

For those distant galaxies, it is possible to use our formalism 
to predict global quantities that can be observed, such as the num¬ 
ber of SNRs expected at a given time and therefore their total X-ray 
luminosity (Rimoldi et ah, in preparation). Exploiting the link be¬ 
tween SNRs and young massive stars, it is also possible to estimate 
the expected SFR in the spheres of influence of quiescent SMBHs. 
These studies, which will be presented in follow-up work, and their 
comparison with observation, can inform theory of nuclear assem¬ 
bly and galaxy formation in general. 
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APPENDIX A: INTEGRALS OE DENSITY IN THE 
EJECTA-DOMINATED STAGE 

Here, we elaborate on the treatment of the integrals over the density 
discussed in Section [3T] The general angle-dependent approach is 
given, as well as the integrations over a sphere for reference. 


AI Angle-dependent integrals 

Beginning with equation we consider explosions either outside 
or inside the break in density. In general, the integrals along R' 
can be partitioned into segments between some radii and 

r>f 

-^upper> 

/ ^upper „ 

R-^R'^ dR', (Al) 

^lower 

which, when summed to compose the full integral from R' — 0 to 
R' = L, give the full mass swept out over the path. The integrals 
are split in such a way to calculate sections that are entirely within 
one of the two possible density gradients lj. 

For explosions outside the break in density, if a shock segment 
has an initial angle tpo > sin“^ (Rh/a), the segment will not cross 
into the region interior to the break, and the integral is fully through 
the u! = 3 region. However, if ipo sin“^ (Rh/a), the element of 
the shock crosses the break. In such cases, the radius R' extending 
from R = a has either one or twcP^ solutions for the intersection 
with the sphere of radius (measured from the SMBH) equal to Rh'. 

RU = 5 (s ± ^s^-d{a^-Rl)^ , (A2) 

where s = 2acos'i/)o. For explosions outside the break, the inte¬ 
gral can be split into three possible regions (with Rh± provided by 
equation lATt : R' < i?b_, or < R' < i?bQ, or R' > i?bQ. 


For^O = sin ^ {Rh/a) and i/iq < sin ^ (Rb/c*). respectively 
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These ranges define the integral limits, where each integral has the 
general form 

r 

Mpart oc a“ / -T- di?'. (A3) 

J {a? + - sR'Y^ 

The same holds for explosions inside the break, except that 
there is only one solution, for the intersection with the surface 
at = _Rb and equation iA3l l has only two sets of limits: R' < 
J?b, and R' > i?b. The general solutions to the angle-dependent 
integrals using these limits are lengthy and are not reproduced here. 


A2 Integrals over a sphere 

Due to the introduction of an axis of symmetry by the offset po¬ 
sition of the sphere, the mass integral may be simplified using 
cylindrical coordinates with this symmetry axis (the cylindrical 
z) aligned on the explosion point. A spherically symmetric den¬ 
sity field (having an origin coincident with that of this cylindri¬ 
cal coordinate system) remains constant over the cylindrical polar 
angle Lp, for a given cylindrical^ a, since the spherical R 

(= x/r2 + 22) is constant. Thus the expression for mass is reduced 
to a double integral. 

In such a coordinate system, for a single power-law density of 
exponent —u), the mass swept up (Ms) by a spherical shock front 
that has expanded through a radial distance L (measured from the 
explosion point, a) 


Ms = 2npaa 


ra + L L^ — (z — a)^ 

7 / 

Ja-L Jo 


r (r^ + At Az. 

(A4) 


For the single density distributions of cu = 1 and 3, this evaluates 
to 


' 2lXpa 


Ms = 


[{a + Lf - \a- L\^ - 6a^L] , 


2TZpa,CL 


In 


Qj L 

a — L 


a 


a > L, 


uj = 1 

uj = 3 . 

(A5) 


In the case a > L for the cj = 1 solution, the term in brackets 
reduces to 2L® such that Mg = 47tpoL®/3, which is the trivial 
a; = 0 solution. 

For a model density in which there is a broken power-law dis¬ 
tribution, the integral over a sphere centred at 2 : = a is less straight¬ 
forward. To avoid introducing complicated integral limits to equa¬ 
tion l lA4t . one approach is to determine the overall quantity by sum¬ 
ming integrals over two density distributions, where the integrand 
for each is restricted using Heaviside step functions, H, that break 
the distribution at specified (spherical) radii, R. Using the cylindri¬ 
cal integral of equation (|M), and taking a density distribution that 
is non-zero between two spherical radii from the origin, R = P 
and R — Q, this effectively constrains the density as 

p^ px [iT +z^-p'^-H (^\/r2 + - q)] . (A6) 

For a two-section broken power-law density with cuin and cuout, the 
total swept-up mass will be Mg = Min + Mout. The solutions for 
u! — 1 and 3 are given below. 


Note that we use R to designate the spherical radial coordinate, and r 
to designate the cylindrical one. 


As an alternative method, solutions forw = 1/2 and 3/2 were 
obtained using integrals of the unconstrained density by instead 
splitting the integrated regions into spherical caps, and adjusting 
the integral limits appropriately. These curves match the behaviour 
of those for lj = 1 and 3, showing agreement between the two 
methods of integration. 


A2.1 Solution for LO = 1 

For an explosion at i? = a and an outer break at i? = c: 

Mi = '^[H{a + L){{a + hf [2 {a + L)- 3a]} 

— Ff (a -I- L — c) {(a -f L — c)^ [2 (a -I- L) -I- c — 3a]} 
-H{\a- L\) {(a - Lf [2 ]a - L] - 3a]} 

+ H {\a — L\ — c) {(]a — L\ — cf [2 ]a — L] -|- c — 3a]}] . 

(A7) 


A2.2 Solution for LU = 3 

In this case, an inner density break at i? = c needs to be applied: 


Ms = 2npoa 
-In 


® H{a + L-c) 

Q, L 


(c — a — L)^ 
2ac 


CL L 


-I- IT (a — L — c) < In 


-f 1 


a — L 


a — L 

+ - 1 


+H{\a-L\-c) 


{c-\a-L\f 

2ac 


■ (A8) 


APPENDIX B: ANALYTIC VOLUME, TIME AND 
VELOCITY EXPRESSIONS IN THE KOMPANEETS 
APPROXIMATION 


Solutions for the volume, time and velocity are given in K92, with 
the given coordinate transformation, for power laws of a; = 0 and 
4. Instead, oj = 1 and 3 are relevant for the presently considered 
density profiles, and we outline here the corresponding solutions. 

In a general radial power-law profile, the solution for the vol¬ 
ume is (see equation[8l as well as equation 16 in K92) 


V = 2n {l-cose)R^AR 

Jr- 



2 (fr 


R^ AR, 
(Bl) 



(iT®r 


— Lu) /2. For cu = 1 and 3, 
w = 1 
w = 3 


(B2) 


are the leading (R+) and trailing (R-) points of the shock (along 
9 — 0, using equation[8}. When comparing the equations for w = 1 
and 3 with the ipo parametrisation, it is important to note that a 
phase shift of 7t in t/o is required to obtain correspondence (of the 
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definition of R- and 7?+, for example) between the two cases. In 
terms of x, this gives the following solutions 


V{x) = Kv X i ^ (1 + 2:2) 

r 


U! = I 
oj = 3 


for 


Kv 


32™® 

3 


(B3) 


(B4) 


Therefore the time can he found from the following integrals of 
equation 


t{x) = Kt X < 


f {1 + u^) dti, tj = 1 


if/ 


(1 + u2) 


(B5) 


dw, Cl) = 3 


for 




256™® po 

3\E (72 - 1 ) ■ 


(B 6 ) 


The u! = 1 result may be written via a power-series expansion 
about a: = 0 , 


ti{x) = Kt 


/2xf^ . X" 


V 5 + 9 



(B7) 


The integral for cu = 3 is divergent for x —^ 1, and is represented 
hy the power series about x = 0 


falx) = Kt 


2 x ®/2 


5 




59x1®/® 


52 



(B8) 


For cj = 1, the value of x = 1 y = t/c signifies the mo¬ 
ment at which the trailing point of the shock reaches the density 
singularity (R- — 0). This is, therefore, also the onset of shock 
self-interactions as other parts of the shock wrap around this point. 
This occurs in a finite time t, and, unlike in the a; = 3 case, so¬ 
lutions still exist heyond x = 1 (there is no hlow-out of the shock 
front, as in the case of cu = 3 where the leading point R+ —>■ oo). 

In terms of x, we therefore have 


t;n(x) = X < 


X® (1 -I- x^) ’ 

(x^ - i f 

, y X® (1 -I- x2) ’ 


U! = 1 

U! = 3 


(B9) 


for 




3XE (72 - 1) 
64™®p(i?) 


(BIO) 


Together, these expressions allow a transformation of the solution 
in terms of x into physical units. We are now in a position to com¬ 
pare the numerical solutions with analytic ones, as is done in Sec¬ 
tion |4^ 


The identity cos (2 cos ^ 2 :) = 2z‘^ — 1 is useful for the solution to 
these integrals. 


APPENDIX C: NUMERICAL TREATMENT OF SHOCK 
SELF-INTERACTIONS 

Here, we outline the treatments of self-interacting segments of the 
shock front in our numerical scheme. In an axisymmetric arrange¬ 
ment, self-interactions can happen in two ways. In the first case, 
any parts of the shock that pass over the axis of symmetry (the x- 
axis in our coordinates) will collide with the complementary part of 
the shock travelling over the axis in the opposite direction. This can 
happen as the shock wraps around the SMBH, where the density is 
at its highest. Other shock self-interactions can he caused hy varia¬ 
tions in the density profile that force flowlines to converge, such as 
near a break in power-law densities (see Fig. lClI for an example). 

We therefore need a routine that can detect self-interactions 
in a general way. A simple implementation would he to examine 
the location histories of the flowlines to determine if any have in¬ 
tersected. In our experience, storing these histories puts too high a 
demand on memory. Instead, two methods for detecting interacting 
regions in the shock front were investigated. 

The first approach is an anticipatory one, which tests the spa¬ 
tial divergence of the velocity of the shock front at each time 
step. This quantity gives an indication of whether portions of the 
shock front are converging. However, care is required in choosing 
the threshold of divergence used to define merging regions of the 
shock, which results in this divergence method being difficult to 
tune. 

In practice, it is more straightforward to use a reactive detec¬ 
tion of self-interactions. The shock front, defined by the positions 
of the flowlines at one point in time, is constrained algorithmically 
to be a simple piecewise linear curve. Any intersections along the 
curve can be detected, as we monitor the ordering of the points. 
Loops arising from intersections of the shock front are removed, 
which is equivalent to an effective merger of all flowlines involved 
in the intersecting loop into a single resultant flowline. 

In all cases of self-interactions, the merged flowlines are re¬ 
placed hy a single flowline with an average of their positions. The 
ejecta mass represented hy the new flowline is taken to he the sum 
of the masses assigned to the previously merging flowlines. 

Fig.lCT] shows an example of regions along a sample shock 
solution in which flowlines are converging and being merged. This 
specific example models fluid elements (flowlines; dashed curves) 
at the shock front (solid lines) colliding due to the change in gradi¬ 
ent of the background density (thick blue line). 
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Figure Cl. Example of flowline intersection and merging due to the pres¬ 
ence of a boundary (thick blue line) between two power-law densities. 
Dashed lines show flowline paths, which lie perpendicular to the shock front 
(sample shock snapshots are given as solid grey lines). In this example, the 
break is between power laws of = 112 and a;out = 3, for an explosion 
around a 10^ Mq SMBH at = 0.6 pc, Ry = 0. 
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